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Breakdown of fluctuation-dissipation relations in granular gases 
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Oh- Abstract 

j_ ■ A numerical molecular dynamics experiment measuring the two-time correlation function of 

the transversal velocity field in the homogeneous cooling state of a granular gas is reported. By 



-C measuring the decay rate and the amplitude of the correlations, the accuracy of the Landau- 



Langevin equation of fluctuating hydrodynamics is checked. The results indicate that although a 
Langevin approach can be valid, the fluctuation-dissipation relation must be modified, since the 



viscosity parameter appearing in it differs from the usual hydrodynamic shear viscosity 
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The hydro dynamic description of ordinary one-component fluids is made by means of 
the local number density, n(r,t), the flow velocity, u(r,t), and the temperature, T(r,t) [l|. 
These fields are determined from some closed macroscopic equations, e.g. the Navier-Stokes 
equations, with the appropriate initial and boundary conditions. On the other hand, the 
most detailed microscopic description of the evolution of the fluid, in terms of the properties 
of the particles forming it, is given by non-equilibrium statistical mechanics. In this latter 
formulation, the hydrodynamic fields appear as the ensemble average of some associated mi- 
croscopic fields, being also possible to describe their stochastic properties, at least formally. 
There is an intermediate level of description, often referred to as the mesoscopic description, 
in which macroscopically fluctuating hydrodynamic fields are considered. 

The first formulation of fluctuating hydrodynamics for equilibrium states was proposed 
by Landau and Lifshitz [1] . Their equations describe the fluctuations of the hydrodynamic 
fields about their equilibrium values, and have the form of linear Langevin equations. The 
noise terms are taken as white and Gaussian with the second moments determined by the 
Navier-Stokes transport coefficients of the fluid. This establishes a relationship between the 
fluctuations of the fields at equilibrium and the dissipation by transport in the fluid, and 
so are known as fluctuation-dissipation (F-D) relations of the second kind [2]. Since the 
formulation of the theory, the question of whether the same idea can be used to describe 
fluctuations in systems that are far away from equilibrium has attracted the attention of 
scientists. For molecular fluids, it seems well established that the F-D relations apply under 
the same conditions as the usual Navier-Stokes hydrodynamic equations 3j. 

Granular fluids are always out of equilibrium due to the kinetic energy dissipation in 
collisions. Evidence has been accumulated during the last years that they provide a unique 
proving ground for non-equilibrium statistical mechanics ideas and techniques. Nevertheless, 
the analogy between granular and ordinary fluids must not be pushed too far. Granular 
materials have many peculiar properties arising from the macroscopic size of the particles 
and the lack of kinetic energy conservation before and after a collision |4J . In particular, the 
validity of the F-D relations can not be assumed a priori. 

Here a definite test of the validity of the F-D relations for granular gases is carried out 
under the most favorable and ideal conditions: the simplest model, state, and hydrodynamic 
field. If the F-D fails in this case, it will also certainly fail for more "realistic" conditions. 
More precisely, the behavior of the transversal component of the flow field in the so-called 



homogeneous cooling state (HCS) of a granular gas is the topic addressed in this paper. 
Conceptually, this problem is similar to the experimental situation investigated in ref. (5|, 
although the granular state considered is different. To start with, the Landau-Langevin 
formalism for the transversal velocity field in an equilibrium ordinary fluid will be briefly 
recalled. 

Consider the local fluctuations of the velocity flow Su(r, t) around the equilibrium value 
u e = 0, and introduce its Fourier representation 5u q (t). This vector can be decomposed 
into the component parallel to q and another component perpendicular to it: Su q (t) = 
Su qt \\q + 6u qt ±(i), with q = q/q and 5u q< ± ■ q = 0. The vector component Su qi ±(t) is the 
vorticity or transversal flow field and obeys the Langevin equation [1] 

d t + (mn e )~ l 7] e q 2 Su qt ±(t) = f q ,±(t), (1) 

where m is the mass of the fluid particles, r\ e the Navier-Stokes shear viscosity of the fluid, 
and f q .±(t) a randomly fluctuating force that is a Gaussian white noise of zero average. It 
verifies the fluctuation-dissipation relation 
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In the above expressions, angular brackets denote average over the noise realizations, the 
index e is used to indicate quantities evaluated at equilibrium, V is the volume of the system, 
I is the unit tensor in the subspace perpendicular to q, and rj e a coefficient that according 
with the Landau theory is the same as the viscosity, r\ e = r] e . This is a main physical 
assertion of the theory. 

There is no equilibrium state for granular fluids. Instead, the simplest state for a freely 
evolving granular fluid is the HCS, mentioned above. At a macroscopic level, it is charac- 
terized by a uniform density n#, a vanishing flow velocity, and a monotonically decreasing 
granular temperature, Tn(t), that obeys the Haff law [6], 

d t T H (t) = -( H (T H )T H (t), (3) 

(h{Th) being the cooling rate accounting for the energy dissipation in collisions. The ques- 
tion addressed here is whether Eqs. (pQ) and ([2]) are also valid to describe fluctuations of 
the vorticity in the HCS of a granular fluid as, in fact, has been assumed sometimes in the 
literature [7|. Of course, if the answer is negative, there is little hope that those equations 
hold for stronger non-equilibrium states, such as vibrofluidized granular media. 



The simplest and most used model for granular fluids is an ensemble of N smooth inelastic 
hard spheres (d = 3) or disks (d = 2) of diameter a [8( . Inelasticity is described by means 
of a constant coefficient of normal restitution a, defined in the interval < a < 1. For this 
model, expressions for the transport coefficients have been obtained by using the inelastic 
Boltzmann-Enskog equation and the Chapman-Enskog procedure. In particular it has been 
found that [3, |lfl| 

Vh = mn H £ v H (t)r]*(a,n H ), (4) 

( H (T H )=v H (t)£- 1 C(a J n H ), (5) 

where £q = (n#<7 d_1 ) _1 is proportional to the mean free path of the particles, vs{t) = 
[2Tfl-(£)/m] ' is a thermal velocity for the HCS, and 77* and (* are dimensionless functions 
of the density and the coefficient of restitution. When Eqs. ([T]) and (J2J) are applied to the 
HCS, the problem arises that the viscosity and the second moment of the noise term are 
time-dependent through the temperature Tfj(t). It is then convenient to introduce dimen- 
sionless length and time scales by £ = r/£ and ds = VH(t)dt/£ , respectively. Moreover, a 
dimensionless velocity field is defined by u(l,s) = 5u(r,t)/vH(t). Then, Eqs. (CQ) and ([2j) 
become 

(d, - C*/2 + V*k 2 ) Wfc, ± (a) = £fc, ± (s), (6) 

(^fc,±(s)^,±(s')>H = V^JV-^Cs - s')<5fc,-fc'^ 2 l, (7) 

with A; = £ q and 1/* = V/^. Moreover, in Eq. (J7]) it is 77* = 77*. For wavevectors k such 
that A_l(&) = C*/2 — ^7*^ 2 < 0, the long time solution of Eq. flS]) reads 

Wh,±(s)=[" d Sl e x ^ s -^£ k , x ( Sl ). (8) 
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Using the above expression and Eq. (JTj), the two-time correlation function of the transverse 
velocity field is easily obtained, 

C k , k ,(s, s') = (u k , ± (s)u> k ,, ± (s')) H = - ^xTk) ^-fe' eAx(fe)(W)| ' ( 9 ) 

valid for s > s' 3> 1. 

In the following, MD simulation results testing Eq. ([9]) will be reported. Consider a square 
cell (d = 2) of side L with periodic boundary conditions, so the minimum value of k allowed 
is k min = 27t£q/L. The value of L has been always chosen such that k min > k_i_ = ((*/2rj*) ' , 
to guarantee that the system be always in the parameter region in which it is expected to 



be stable with respect to vorticity fluctuations, i.e. \j_(k m in) < and Eq. ([9]) holds for all 
the allowed values of k. In the simulations, fluctuations corresponding to k = k m i n e x , where 
e x is the unit vector in the direction of one of the sides of the cell, have been measured. 
Therefore, 0Jh,± — ^k,L^y, with e y being an unit vector perpendicular to e x . The particle 
representation of the transversal velocity field considered here is 

N 

^ m ^M = E^ 2mXm/L v y M (io) 

where the components of the position and velocity of particle i at time t have been denoted by 
Xi(t),Yi(t) and V X) i(t), V V) i(t), respectively. In the MD simulations an event driven algorithm 
and the steady state representation of the HCS |lll.ll2| have been employed. The fluctuations 
of the transverse velocity turned out to be always Gaussian within the statistical errors. 
Moreover, the imaginary part of the two-time correlation function was found to vanish, as 
expected by symmetry considerations and in agreement with Eq. (Q. The restriction to 
the parameter region in which the HCS is stable implies that the size of the system and, 
therefore, the number of particles at a given density, is bounded from above. Consequently, 
relevant size effects show up at low densities. In Fig. [TJ 
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is plotted as a function of s — s' for a = 0.8, n#cr 2 = 0.01, and several values of N, as 
indicated. In all cases, an exponential decay is observed, in agreement with Eq. (jUJ), that 
predicts $(s, s') = exp [Aj_(fc m i n )(s — s')]. From the data, values for \±(k min ) are obtained. 
Of course, similar results are found when the labels of the two axis are interchanged. The 
simulations also provide the value of the cooling rate, that is independent of iV and in good 
agreement with the theoretical prediction given in [9|, [lOj. Then, by subtracting C*/2 from 
the obtained values of X±(k m i n ), simulation results for the parameter rf follow. From the 
analysis of the data, it follows that the value obtained depends on the number of particles 
iV (size L) employed in the simulation. The theoretical prediction corresponds to the formal 
limit iV — > oo, L — ► oo, with n^ constant. In Fig. [21 the values of if obtained from the 
data in Fig. [T]are plotted as a function of iV^ 1 . A linear fit leads to the result liniTv^oo if = 
(1.059 ± 0.004)?7*, where rf e = (2^/2n)~ 1 is the low density (Boltzmann) elastic limit. The 
theoretical prediction for the Navier-Stokes shear viscosity obtained from the Boltzmann 
equation in the first Sonine approximation is vf = 1.057??*. Moreover, if a "size-dependent 
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FIG. 1: MD simulation results showing the exponential decay of the two-time correlation function 
of the transversal velocity for a = 0.8 and n#<r 2 = 0.01. 

viscosity" is measured by considering the decay of a small macroscopic sine perturbation 
of the velocity field [13| in systems of different sizes, the values are in good agreement 
with those reported in Fig. [2j Therefore, it can be concluded that, for the system being 
investigated, Eq. (J5J) predicts accurately the decay of the correlation function, with rj* being 
the (inelastic) Navier-Stokes shear viscosity. 

Focus next on the amplitude of the fluctuations. From Eq. (Q, it follows that 
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V* 2 n*k 2 
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(12) 



Since X±(k m i n ) was determined from the time decay of the correlation function, the coefficient 
rf can be computed from the measurements of the above quantity. Again, a clear dependence 
on the number of particles iV is observed in the results reported in Fig. [31 that correspond to 
the same system as in Figs. [T]and[2j By means of a linear fit, the asymptotic N —>■ oo value 
rf = 1.29 ± 0.02 is found. It differs from both the shear viscosity of the system measured 
by the procedures described above and the theoretical prediction in more than 20%. 

The above analysis has been repeated for different values of the coefficient of restitu- 
tion, keeping the density small to allow comparison with the kinetic theory predictions. In 
particular, the results plotted in Fig. @] correspond to the case previously discussed and to 
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FIG. 2: Values of the parameter rf /rj* = 2\/2in]*, obtained from the exponential decay of the 
correlation function of the transversal velocity field, as a function of the inverse of the number of 
particles N. The density and coefficient of restitution are the same as in Fig. [TJ Symbols are from 
MD simulations and the solid line is the best linear fit. 

a = 0.85, 0.90, and 0.95, with n H a 2 = 6 x 10~ 3 , 0.018, and 0.02, respectively. It is observed 
that the discrepancy between the shear viscosity rj and the coefficient rj characterizing the 
second moment of the fluctuations, increases as the inelasticity increases, vanishing in the 
elastic limit, as required. In the simulations, attention has been paid to assure that k m i n 
is neither too large nor too small. In the former case, the Navier-Stokes approximation is 
expected to fail, while in the latter instability effects show up. 

Two main physical conclusions follow from the study presented here. The Landau- 
Langevin fluctuating hydrodynamic equations cannot be directly extended to granular fluids. 
But, on the other hand, it seems possible to separate slow and fast degrees of freedom and 
to model granular systems in the spirit of the Langevin approach. A similar conclusion was 
reached from the study of the fluctuations of the total energy of the system [14J]. The open 
and relevant question is how to express the amplitude of the noise source (i.e. rf) in terms 
of the macroscopic properties of the system. Answering this question requieres to derive 
Langevin-like equations for the hydrodynamic field of granular fluids stating from kinetic 





1 

o 


i ' 


i ' 


i 




1.15 


>. < 


> 






- 


rf/r|* e 










- 


1.1 










- 









N. O 




- 


1.05 










- 


o!c 


1 


1 


1 


1 




)004 


0.0006 


0.0008 


0.001 





1/N 



FIG. 3: Values of the parameter rj* /rj* = 2y/2irrj* , obtained from the amplitude of the fluctuations 
of the transversal velocity field, as a function of the inverse of the number of particles. The density 
and coefficient of restitution are the same as in Fig. [TJ Symbols are from MD simulations and the 
solid line is the best linear fit. 

theory or statistical mechanics descriptions. 
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